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ABSTRACT 

We introduce a fast Markov Chain Monte Carlo (MCMC) exploration of the astrophysical parameter 
space using a modified version of the publicly available code CIGALE (Code Investigating GALaxy 
emission). The original CIGALE builds a grid of theoretical Spectral Energy Distribution (SED) 
models and fits to photometric fluxes from Ultraviolet (UV) to Infrared (IR) to put contraints on 
parameters related to both formation and evolution of galaxies. Such a grid-based method can lead 
to a long and challenging parameter extraction since the computation time increases exponentially 
with the number of parameters considered and results can be dependent on the density of sampling 
points, which must be chosen in advance for each parameter. Markov Chain Monte Carlo methods, 
on the other hand, scale approximately linearly with the number of parameters, allowing a faster 
and more accurate exploration of the parameter space by using a smaller number of efficiently chosen 
samples. We test our MCMC version of the code CIGALE (called CIGALEMC) with simulated data. 
After checking the ability of the code to retrieve the input parameters used to build the mock sample, 
we fit theoretical SEDs to real data from the well known and studied SINGS sample. We discuss 
constraints on the parameters and show the advantages of our MCMC sampling method in terms of 
accuracy of the results and optimization of CPU time. 

Subject headings: galaxies: fundamental parameters - methods: data analysis 



1. INTRODUCTION 

The spectral energy distribution (SED) of galaxies de- 
pends on many physical processes related to the emis- 
sion from different stellar populations, absorption and 
re-emission from dust and gas and possible presence of 
Active Galactic Nuclei (AGN). Each process has been 
studied by many authors; libraries of stellar population 
models (Fioc & Rocca-Volmerange (1997), Bruzual & 
Chariot (2003), Maraston (2005)), fitting curves for dust 
emission (Calzetti et al. (1994, 2000), Witt & Gordon 
(2000)), studies of emission of dust grains (Chary & El- 
baz (2001), Dale & Helou (2002), Lagache et al. (2003, 
2004), and Siebenmorgen & Kriigel (2007), Silva et al. 
(1998), Dopita et al. (2005), da Cunha et al. (2008)) 
are the basis of sophisticated fitting codes which derive 
physical parameters such as stellar mass, star formation 
rate, dust luminosity and so on. 

Many parameters are usually necessary to describe these 
processes and model theoretical SEDs of galaxies. A grid 
of theoretical SED models is usually built and fitted to 
the data and statistical properties are derived for the pa- 
rameters of interest. A big drawback of any grid-based 
method is that, for any fitting process, the time to build 
models grows linearly with the number of models and 
then about exponentially with the number of parame- 
ters involved: such approaches are difficult to implement 
for complex models involving a sufficiently large num- 
ber of parameters or when a fine sampling of the pa- 
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rameter space is necessary in order to retrieve statis- 
tically robust results. In the past few years, Markov 
Chains Monte Carlo (MCMC) techniques have started 
being widely used in physics. In cosmology, parameter 
estimation from cosmic microwave background data with 
MCMC methods has been introduced in Christensen 
et al. (2001) and has been implemented in the pub- 
licly available code cosmomc (Cosmological Monte Carlo, 
Lewis & Bridle (2002) j3; in astrophysics, an MCMC ap- 
proach to the stellar population syntesis modeling has 
been introduced in Conroy et al. (2009). 
Here we use cosmomc as a generic sampler and we inter- 
face it to the publicly available code CIGALE (Code 
Investigation GALaxy Emission, Noll et al. (2009)) in 
order to allow a fast and accurate evaluation of the mul- 
tidimensional parameter space probed by this code 0. 
The main advantage of this method is that the comput- 
ing time to fit the data scales approximately linearly with 
the number of parameters involved, allowing the user to 
consider complex models with many parameters for only 
small additional computational time. MCMC techniques 
allow to probe also the shape of the probability distribu- 
tion, giving far more information than just best fit and 
marginalized values for the parameters. 
The paper is organized as follows; in the next section we 
briefly describe CIGALE, introducing the main parame- 
ters used in the subsequent sections. We then explain the 
MCMC technique implemented in the modified version 
of CIGALE, which we caU CIGALEMC. We test our code 

http : // cosmologist . inf o/cosmomc/ 
^ http://www.oamp.fr/cigale/ 

^ During the completion of this work we noticed that Acquaviva 
et al. (2011) have performed a similar work in the context of the 
code GALAXEV developed by Bruzual & Chariot (2003). 
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using a mock sample of 62 galaxies already used in Gio- 
vannoli et at. (2011) and we apply it to a real galaxy sam- 
ple with data from the Spitzer Infrared Nearby Galaxy 
Survey (SINGS, see Kennicutt et al. (2003)). We al- 
ways consider a flat cosmological model with fi™ = 0.3, 
Qa = 0.7 and Hq = 70kms"^Mpc"\ Finally we give 
our results and conclusions. 

2. THE CODE CIGALE 

CIGALE calculates a grid of theoretical SEDs and fits 
to observational input data constitued by photometric 
filter fluxes ranging from UV to IR. For a detailed de- 
scription of the code and its application to real data, we 
refer the interested reader to these papers (Burgarella et 
al. (2005), NoU et al. (2009), Giovannoh et al. (2011), 
Buat et al. (2011)). In the following, we briefly sum- 
marize its main characteristics and the basic parameters 
used in the next sections. Our notation follows the one 
introduced in Giovannoli et al. (2011). 

2.1. Stellar populations and star formation rate 

CIGALE combines both old and young stellar popula- 
tions using single stellar populations of Maraston et al. 
(2005) or PEGASE (Fioc & Rocca-Volmerange (1997)). 
In this paper we will only use Maraston models; we as- 
sume star formation histories (SFH) with either exponen- 
tially decreasing star formation rate (SFR) in function of 
time t ("r models"), as: 

SFRoid(t) = SFRo^oid • (1) 

(*-*2) 

SFRyoung(t) = SFRq young ' 6 ^2 

or "box models" characterized by constant SFR over a 
limited period of time; in this case the instantaneous 
SFR at look-back time t' = is given by the galaxy 
mass divided by the age t of the population, i.e. Mgai/t. 
Labels 1 and 2 refers to the old and young stellar pop- 
ulations while Ti_2 and ti^2 (both is units of Gyr) are 
their e-folding time and age respectively. The two stellar 
populations are linked and weighted through their mass 
fraction; the parameter /ysp represents the fraction of 
the young stellar mass over the total mass, so that the 
total instantaneous SFR (output parameter of CIGALE) 
is expressed as: 

SFR - (1 - fySp)SFRoid(t) + fySpSFRyou„g(t). (2) 

2.2. Absorption and emission by dust and gas 

In CIGALE, the absorption of star light by dust is 
described by a Calzetti attenuation curve (Calzetti et 
al. (1994, 2000)); possible modifications of the curve 
include both the addition of a UV bump at about 
2175A(Fitzpatrick & Massa (1990; 2007); NoU et al. 
(2009)) and the change of the slope through the mul- 
tiplication by a power law (A/Ay)*, where Ay = 5500A 
is the reference wavelength for the V filter; both the am- 
plitude of the bump and the slope S are free parame- 
ters of the code. The attenuation correction is applied 
to both stellar populations individually using the visual 
attenuation parameter of the young stellar populations 
( AySp , in units of magnitudes) and a reduction factor of 
the attenuation for the old model (/y) as free parame- 



terfl The code also calculates 3 additional parameters, 
derived from the final model SEDs; Apuv and Ay are 
defined as the effective attenuation factors in magnitudes 
at 1500 ± lOOA and 5500 ± lOOA while the age tD4ooo is 
derived from the D4000 break (see Balogh et al. (1999)) 
of the unreddened SED for a single stellar populatior{3. 
Dust emission in the IR is taken into account using 64 
templates of Dale & Helou (2002). These models are 
parametrized by a, the power law slope of the dust mass 
over heating intensity, defined as follows: 

dMd{U) = U-'^dU, (3) 

where Md{U) is the dust mass heated by a radiation field 
of intensity U. 

Bolometric and dust luminosities (Lboi and Ldust respec- 
tively) are derived from all the basic parameters (see Noll 
et al. (2009)) and dust emission due to non-thermal 
sources such as AGN can also be added; the fraction 
/agn of dust luminosity Ldust (in Lq) due to an AGN is 
estimated using AGN templates from Siebenmorgen et 
al. (2004). 

The spectral line correction due to interstellar gas is per- 
formed as in Noll et al. (2009): for the optical band, 
empirical line templates are taken from the Kinney et al. 
(1996) starburst spectra while for the UV we use tem- 
plates derived from SEDs presented in Noll et al. (2004). 
A correction for the redshift-dependent absorption of the 
intergalactic medium shortward of the Lya line is also in- 
cluded using the algorithm of Meiksin (2006). 

2.3. Comparison with data 

A grid of theoretical photometric fiuxes is calculated 
at the redshift of the objects considered and a Bayesian 
analysis is performed through the calculation of the 
of each model: 

x'(Mgai) = ^^"^'"^'""t ' ~ ^""'''^^ = -21n(L); (4) 

i—\ ^obs,z 

here the galaxy mass Mgai (in Mq) is treated as a free 
parameter, /mod,i and /obs.i are the theoretical and ex- 
perimental fluxes respectively, the statistical photometry 
errors are considered in the term (Tobs,i and L is the nor- 
malized likelihood function. 

3. MCMC TECHNIQUE AND COSMOMC 

In Bayesian inference, the posterior probability of the 
parameters {9) of a model in the light of the observed 
data (d) is given by: 

P(.|J) = »5^; (5) 
P{d) 

here P(d|6') = L(6') is the likelihood of the data given 

the model, P(0) is the prior on the parameters, which 
quantifies our a priori knowledge of the parameters and 
P(d) (called Evidence) is a normalization factor. In 

^ The parameter fy has been labelled /att in some previous 
papers. 

* In its current version, CIGALE directly outputs the dust-free 
D4000 break, see Buat et al. (2011). 
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Fig. 1. — Samples from the posterior distribution for a test 
galaxy; the high density of points in the parameter space corre- 
sponds to large values of the posterior. Units are MQyr"-' , M©, Lq 
for SFR, Mstar, Ldust respectively. 

our case, d represents the SED of each galaxy while 6 
represents the astrophysical parameters of CIGALE, as 
{6}i = {Ti,t2, fySP, ■■■}■ An MCMC sampler provides 
an efficient way to explore the posterior distribution and 
ensures that the number density of samples is asymptot- 
ically proportional to the probability density. 

3.1. Metropolis-Hastings algorithm 

The code cosmomc uses the Metropolis-Hastings algo- 
rithm to generate samples; each chain moves according to 
a transition probability T(9i,9i^i) which is determined 
so that the Markov Chain has a stationary asymptotic 
distribution equal to the posterior distribution P{9) that 
we want to sample from. Given an arbitrary proposal 
density distribution q{9i,di-^-i) to propose a new point 

9i+i when the chain is at the point 9i, the probability of 
transition is given by /3: 



iin|l, 



Pi9,)qi0^,0^+l) 

SO that 

T{9,, 9,+,) = mJ,+,)q{9,, 9,+,). (7) 
This ensures that the detailed balance holds: 

P{9,+i)T{9,+i,9,) = P0i)T0„ 9,+i) (8) 

and that the distribution converges to P(6'). In practice, 
a random number a; S [0 : 1] is generated in the process 
of moving from 9i to so that the new point is 
accepted if /3 > x. This ensures that each point of the 
chain depends only on its predecessor; in this sense the 
chain is a Monte-Carlo Markov process. 

3.2. Comparison with grid-based methods, burn in and 
convergence diagnostics 

As an illustration of the sampling mechanism, in Fig- 
ure 1 we plot samples from the posterior distribution 



for a MCMC run with a test galaxy taken from a mock 
sample at redshift z ^ 0.7 (see the following section for 
details); the number density of samples in the plane is 
proportional to the probability density of these two pa- 
rameters. The dust luminosity Ljust strongly depends 
on the SFR and the two parameters are degenerate, as 
shown by the colours in the figure. This plot clearly 
shows the efficiency of this MCMC method. In the grid- 
based approach, the parameter space is sampled in the 
same "blind" way for high and low values of the poste- 
rior: this can be an issue for both reliability of results 
and computation time, as also pointed out in Noll et al. 
(2009) and Acquaviva et al. (2011). In fact, local min- 
ima and degeneracies between parameters can be easily 
missed or undersampled without a good a priori knowl- 
edge of the parameter space; the oversampling of an ill- 
constrained parameter can also lead to a slight degrada- 
tion of the estimates of well constrained parameters and 
many points can be generated in a region where the pos- 
terior is low, resulting in a waste of CPU time. This is not 
the case when MCMC chains are used because each chain 
"learns" where to move in the parameter space through 
the Metropolis-Hastings algorithm so that the density 
of samples is proportional to the posterior distribution. 
Degeneracies between parameters are more easily found, 
especially if many chains, starting from different regions 
in the parameter space, arc used. In other words, the 
CIGALEMC user needs to specify the prior parameter space 
(number of parameters and their limits) but not the den- 
sity of points for each parameter. In the following section 
we will provide a comparison of CPU time between the 
original CIGALE and CIGALEMC when evaluating physi- 
cal properties of a mock sample. The code also calculates 
the covariancc between various parameters so that an ini- 
tial run can be made and the covariance matrix obtained 
can be used to improve the efficiency of sampling for sub- 
sequent runs. 

Since each MCMC chain starts at a random position in 
the parameter space, it will take a little time before the 
chain equilibrates and starts sampling the posterior dis- 
tribution. This period of initial convergence is called 
burn in period and the first burn in points of each chains 
will be discarded when doing any statistical analysis. In 
order to obtain uncorrelated samples of the posterior 
each chain is also "thinned" by using only occasional 
points of it; the thinning factor varies according to the 
number of parameters involved and it is tipically in the 
range 25-50. The code allows to choose the burn in frac- 
tion of the chain we want to discard and automatically 
thins out the chains. 

In our analyses we won't use any initial covariance matrix 
for the parameters and, in order to be conservative, we 
compute statistical quantities using only the second half 
of each chain. In Figure 2 we plot the points of a MCMC 
chain for a galaxy sample in the plane Ldust vs fAGN ; the 
chain reaches the sensitive region of the parameter space 
after only a few "burn in" points characterized by very 
low values of Ldust- Having a set of samples from the 
full posterior distribution, it is possible to calculate sta- 
tistical quantities for the parameters of interest. Since 
the number density of sampling points is proportional to 
the posterior density, it is easy to calculate mean values 
and marginalized 1-dimensional distribution for each pa- 
rameter 9^ by simply counting the number N of samples 
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Fig. 2. — A Monte Carlo Markov chain in the 2-dimensional 
parameter space Ljuat vs Jagn- The chain starts in a region where 
the likelihood is low ("burn in" points with L(just ~ 8.8L0) and 
quickly reaches the most sensitive region in the parameter space. 



within binned ranges of parameter values: 



1 ^ 
-T 

N ^ 



(9) 



while this is much more difficult in the context of the 
numerical grid integration because the calculation time 
grows exponentially with the number of dimensions. 
In order to be sure that MCMC chains are efficiently sam- 
pling the posterior distribution (and then obtain robust 
statistics for each parameter) it is important to check 
their convergence. The code cosmomc provides two con- 
vergence criteria for runs with one single chain (Raftery 
& Lewis, 1992) and with multiple chains (Gclman & Ru- 
bin, 1992). In the following analysis we will run multiple 
chains using the Gelman & Rubin diagnostic which is 
characterized by the "variance of chain means" /"mean 
of chain variances" parameter R; |R — 1| < 0.03 is usu- 
ally enough to reach convergence and stop the chains. 
In this work we use cosmomc as a generic MCMC sam- 
pler and we link it to CIGALE in order to allow a faster 
exploration of the astrophysical parameter space. Our 
modified CIGALEMC code will be publicly available very 
soor0 and, since it is based on cosmomc for sampling op- 
tions, convergence criteria and statistical quantities pro- 
vided, we refer the reader to the websit43 and to Lewis 
& Bridle (2002) and references therein for a detailed ex- 
planation of the code and MCMC methods in general. 

4. ANALYSIS OF A MOCK SAMPLE 

We test our CIGALEMC code with a mock sample already 
used in Giovannoli et al. (2011). We consider 62 arti- 
ficial galaxy SEDs corresponding to Luminous InfraRed 
Galaxies (LIRGs) at redshift z ~ 0.7 and obtained by 
varying the input parameters of CIGALE in the follow- 
ing ranges: [0.1 Gyr < n < 10 Gyr, 0.025 Gyr < t2 < 
0.7Gyr,0 < fysp < 1,0.6 mag < Aysp < 2.1 mag, < 
fv < 1,1.0 < a < 2.5,0 < Jagn < 0.3]; such a wide 
multidimensional parameter space allows to model very 
different spectral energy distributions of galaxies charac- 

® http://www.Ocmip.fr/cigale/. 

http : // cosmologist . inf o/cosmomc/ 



terized by a wide range of possible star formation histo- 
ries, absorption and emission by gas and dust and AGN 
contamination. All galaxies are based on a Salpeter ini- 
tial mass function; the age of the old stellar population 
is fixed at 7 Gyr, we consider a constant star formation 
rate for the young population model {t2 — 20 Gyr) and 
we do not add any modification to the original Calzetti 
attenuation curve (no UV bump and 6 = 0). 
Theoretical fluxes are calculated in the following 17 
bands from UV to IR: 0.231/im for GALEX, [0.35-0.36- 
0.46 - 0.54 - 0.65 - 0.87 - 0.90 - 1.2 - 1.6 - 2.1]^m (cor- 
responding to MUSYC bands), [3.6 - 4.5 - 5.8 - 8.0]jum 
(corresponding to IRAC photometry) and [24 — 70] /im 
(corresponding to MIPS photometry). We add a gaus- 
sian distributed uncertainty a to each theoretical flux; 
its value is 10% of the corresponding flux, a reasonable 
choice for the measurements considered above. For each 
galaxy we run 8 chains with initial positions randomly 
chosen in the parameter space determined by the follow- 
ing set of input parameters: 



{ri, t2, fySP, AySP, fv, a, fAGN, Mgal}. 



(10) 



from which statistical quantities of interest can be calcu- 
lated for the following set of derived parameters : 



{SFR, tD4000i Lbolj Ldust, AFUV,Av,Mstar}. 



(11) 



We assume flat priors in the following parameters 
space: [0.1 Gyr < n < 10 Gyr, 0.025 Gyr < t2 < 
2Gyr,0 < fySP < 1,0 mag < AySP < 3mag,0 < fv < 
1,0.5 < a < 3,0 < fAGN < 1,8M0 < Mg^i < UMq]; 
chains are stopped when the Gelman & Rubin R-1 pa- 
rameter is |R — 1| ~ 0.03. 

First of all, we want to check that our results are statis- 
tically in agreement with the input values for the mock 
sample. As a tool derived from cosmomc, CIGALEMC 
allows to calculate and plot the mean likelihood and 
marginalized distribution for each parameter. 

The marginalized distribution in a given direction of 
the parameter space d = h{9) (where h{d) is the pro- 
jector operator in one of the parameters considered, as 
h{6) = 6i) is proportional to the number of samples at d 
and it can be expressed as: 



P{v) = MiP,v) = / d9P{e)Sih{d) - v), 



(12) 



where P{0) is the posterior distribution. Assuming flat 
priors on 9, the mean likelihood of samples with d = h(6') 
can be expressed as: 



{P{v:h{9) = v)) 



jdep{e)^6{h{e)-v) _ m{p^,v) 



J de p{e)d{h{e) - v) m (p, v) ' 

(13) 

If P{0) is a multivariate Gaussian distribution it is 
possible to demonstrate that both mean likelihood and 
marginalized distribution are Gaussian and proportional 
so that they look the same: differences in these distri- 
butions will be then a signal of non-Gaussianity which 
can be either intrinsic or due to parameters not very well 
constrained. In Figure 3 we show, for an artiflcial galaxy 
of the mock sample considered, both marginalized distri- 
butions (black solid lines) and mean likelihood (dotted 
lines) for some parameters of interest: we can see that 
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Fig. 3. — Mean likelihoods, computed over the whole chains, for binned parameter values (dotted lines) and marginalized distributions 
(black lines) of some parameters of interest for an artificial galaxy of the mock sample considered; the blue vertical lines show the input 
parameter values used. Poorly constrained parameters, such as ri, i2, fySPt fv ^nd a are clearly visible. Luminosities and masses are in 
units of solar luminosities and solar masses respectively. 
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Fig. 4. — Distribution of the variable X for some parameters considered; X values are compatible with at 68% c.l.. 
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Fig. 5. — Two-dimensional marginalized distribution showing the 
68% and %95 c.l. contours for fySP and t2 for a galaxy sample; 
the parameter t2 is unconstrained by the data and it is partially 
degenerate with fySP 



Til ^2, fySP fv a-nd a arc not very well constrained by 
the code. Similar results have been found in the analyses 
by NoU et al. (2009), Giovannoli et al. (2011) and Buat 
et al. (2011). 

In order to study the goodness of the fit in a quantita- 
tive way for the whole sample of galaxies we introduce 
the quantity: 

x'-^E^; (14) 



here i runs for all the parameters considered by the code 
(so that {X}' = {Ti,t2,fysp, ...}), j runs for the iV 62 
objects of the mock, O' is the set of the best fit values ob- 
tained as output from CIGALEMC, P is the vector of input 
parameters and cr' is the vector of the 68% c.l. marginal- 
ized uncertainties. As we can see from Figure 4, all X 
values are compatible with X = 0, which means that the 
code is able to find the best fit values of the parameters 
with high confidence. However, some distributions are 
slightly skewed: this can be due to chains being stuck in 
local minima of the likelihood function so that, in some 
cases, the best fit found does not correspond to the input 
value. This situation is typical when parameters are not 
well constrained (as for ^2) or in the presence of strong 
degeneracies between parameters. As an example, Fig- 
ure 5 shows the 68% and 95% c.l. in the plane fySP vs ^2 
for one sample galaxy of the mock: as we can see, limits 
on these parameters are not very strong; we also notice 
a partial degeneracy for high values of fySP which is due 
to the fact that these two parameters affect the galactic 
SFR in the same way since we consider constant SFR for 
the young model: 



SFR ~ fySP 



M, 



(15) 



In order to check the reliability of our MCMC algo- 
rithm we tested the parameter uncertainty estimation. 
We created and fit the SEDs of 100 artificial galaxies in 
21 bands from IR to UV, built with the same input pa- 
rameter set and with a 10% scatter in flux, corresponding 
to the photometric error considered. We checked that in 



Para, met ers 


mock sample 




0.33 


'-'-'£) iU -'•-'^SLar 


0.93 


loP'in Tj\-. 1 


0.92 




0.88 


logio SFR 


0.81 




0.26 


fySP 


0.30 


AySP 


0.81 


fv 


0.70 


a 


0.65 


Jagn 


0.71 



TABLE 1 

Estimation of the linear correlation coefficient of 
Pearson between the exact value and the value estimated 
by the code cigalemc for some parameters of the mock 
catalogue. 



general we are able to find the true values within the re- 
gion allowed at 68 (95)% confidence 68 (95)% of the time 
within the Poisson fluctuation error for different choices 
for both the input parameters and the prior ranges al- 
lowed in the fit. 

Finally, we checked the consistency of our results by cal- 
culating the Pearson correlation coefficient r between the 
input values of the parameters used to generate the mock 
and the best fit values found with CIGALEMC. The Pear- 
son correlation coefficient quantifies the amount of cor- 
relation between two variables X and Y and it assumes 
values in the range r G [—1 : 1] (see Cohen (1988)). For 
samples Xi and Yi, it can be written as: 



TJ^=,{x,~x){y,~y) 



(16) 



here X and Y denote the mean values of Xi and Yi. The 
interpretation of the strength of the correlation can de- 
pend on the context; however, a widely used standard 
introduced by Cohen (1988) considers Pearson values 
\r\ > 0.5 as a signal of large correlation between vari- 
ables, while \r\ < 0.3 denotes poor correlation. Again, 
as we can see from Table 1, some parameters (ri, ^2, 
fySp) have small values of r; this is expected since these 
parameters are mostly unconstrained; however, even if 
the exact values of the input parameters are not found, 
they always fall inside the statistically significant region 
of the parameter space. We also explicitely checked that 
results do not significantly change when increasing the 
prior space for some parameters not very well constrained 
as Ti, t2 and a. Finally, we tested, for a galaxy of the 
mock, how results change when considering a gaussian 
prior, with different central values and width a — 0.2, on 
one of the most unconstrained parameters, ^2! we found 
that results on the other parameters are always statisti- 
cally consistent respect to the choice of a flat prior on t2- 
In general, it is a good idea to choose the largest possi- 
ble prior space for unconstrained parameters, in order to 
avoid possible biases due to the prior choice. 
It is useful to compare the performance of CIGALE and 
CIGALEMC in terms of computing time, especially since 
computation can become prohibitive for any grid-based 
method if the number of parameters involved is suffi- 
ciently high. The CPU time required to obtain conver- 
gence of the chains for each galaxy mainly depends on 
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both the quahty of the data and the number of param- 
eters considered. Running 8 chains in parallel (each one 
on a 2.40 GHz Intel Xeon E5530), for each galaxy of 
the mock we typically reach a good convergence with 
^ 20000 points for each chain, which means a total of ^ 
160000 points. The grid buiU in Giovannoh et al. (2011) 
to analyze the same sample with the same number of 
free parameters and bounds contained ~ 3.5 • 10^ points; 
this means a gain of order ~20 in efficiency in the esti- 
mation of the parameters but a more dramatic efficiency 
can be easily reached when we need to use either more 
parameters or a fine sampling in a given direction of the 
parameter space or both. The average CPU time to reach 
a good convergence for a galaxy of this mock is about 35 
seconds, which translates in 280s of total CPU time. 

5. ANALYSIS OF REAL DATA: THE SINGS SAMPLE 

We now use CIGALEMC to infer physical properties of 
the well known SINGS (Spitzer Infrared Nearby Galaxy 
Survey; Kennicutt et al. (2003)) sample. In order to 
make a comparison with results obtained using the grid- 
based CIGALE, we use the same 39 galaxies already con- 
sidered in Noll et al. (2009) with the same spectral cov- 
erage: GALEX FUV (- 1500A) and NUV (- 2300A) 
filters (Gil de Paz et al. (2007)), 2MASS data for J,H, 
K, (Jarrett et al. (2003)), IRAC and MIPS fihers for 
[3.6, 4.5, 5.8, 8.0, 24, 70, 160]/im (Dale et al. (2005)), Dale 
et al. (2007, 2008) optical data for B, V, R and I bands 
corrected as in Muhoz-Mateos et al (2009) and fluxes 
from u\ g' , r', i' and z' filters of SDSS (Stoughton et al. 
(2002)); Dale et al. (2007, 2008) optical data are only 
used for 14 galaxies for which SDSS data are not avail- 
able. The mean photometric relative uncertainties for 
the bands considered are shown in Table 2; in particular, 
very small uncertainties are associated with the 2MASS 
bands. We performed a preliminar run, realizing that 
the hard constraints coming from this filter set did not 
allow the code to properly fit for the other filters. Since 
these uncertainties do not take into account for possible 
calibration errors and since a systematic offset can also 
affect the CIGALE theoretical models, we decided to be 
conservative and, following Noll et al. (2009), we add a 
5% uncertainties in quadrature for each filter set. 

In our analysis we assume the following ange of vari- 



Filters 


Rel. errors 


GALEX FUV, NUV 


15% 


Dale ct al. B, V, R, I 


16% 


SDSS u , g , r , i, z 


3% 


2MASS J, H, Ks 


1% 


IRAC 3.6, 4.5, 5.8, 8.0 fim 


11% 


MIPS 24 /J m 


5% 


MIPS 70 fj, 


7% 


MIPS 160 fim 


13% 



TABLE 2 

Mean photometric uncertainties for our SINGS sample. 



ation for a set of 9 astrophysical parameters: [0.1 Gyr < 
Ti < 10 Gyr, 0.025 Gyr < t2 < 2Gyr,0 < fysp < 
1,-0.5 < 5 < 0.5,0 mag < Aysp < 5mag,0 < fv < 
1,0.5 < a < 3,0 < fAGN < l,8M0,Mgai < 13Mo]. 
We keep fixed both the age of the old stellar population 



{ti — 10 Gyr) and the e-folding time for the young stellar 
population (t2 = 20 Gyr): these parameters are not well 
constrained by the data so that fixing them does not alter 
the fit. Finally, We only consider models with Salpeter 
initial mass function and solar metallicity; metallicity 
measurements are quite uncertain due to many limit- 
ing factors (Noll et al. (2009), Moustakas & Kennicutt 
(2006) and references therein); Noll et al. (2009) checked 
the influence of different assumptions for the metallicity, 
concluding that deviations in the properties of the galax- 
ies are within the uncertainties, which tend to increase 
by to 20% when half or double of the solar metallicity 
are considered. The only exception is the absolute value 
of t_D4000j because of the well-known age- metallicity de- 
generacy (e.g., Kodama & Arimoto, (1997)). 
Our reference AGN model has L ~ IO^^Lq, R ~ 125pc 
and Ay = 32 for the luminosity of the non-thermal 
source, the outer radius of a spherical dust cloud cov- 
ering the AGN and the amount of attenuation in the 
optical caused by the cloud respectively. 
Our findings can be summarized as follows: 

• Very good constraints arc derived for the AGN 
fraction of all sources {Jagn < 0.10) except for 
NGC0584 and NGC1404 for which fAGN < 0.6 
at 68% C.I.. We note that the flux at 160^m 
for NGC0584 can be contaminated by some fore- 
ground/background emission and the most recent 
Herschel data conclusively confirm that a back- 
ground source contaminates both fluxes at 70 and 
160/im for NGC1404 (Daniel Dale, private commu- 
nication); the low quality of these data is respon- 
sible for the big uncertainties obtained for other 
parameters as SFR Lboi Ldust Apuv ! in Figure 6 we 
plot the 2-dimensional marginalized distribution 
for SFR vs Jagn for NGC1404: the double-peaked 
likelihood function is most probably an artifact due 
to the low quality of data for this source. In gen- 
eral, we note that that NGC0584 and NGC1404 are 
elliptical galaxies which tend to have very weak but 
warm dust emission; it is not surprising that they 
show an apparently high AGN fraction. We de- 
cided not to use these sources in the rest of our 
analysis. 

• We are not able to put strong constraints on "phe- 
nomcnological" parameters as ri, t2 and 5; limits 
on these parameters depend on the assumed prior 
range. The poor determination of 5 is essentially 
due to the low number of data in UV. In general, 
both 5 and the possible UV bump allowed by the 
code are very difficult to constrain with few only 
broad band data (Buat et al. (2011) and Buat 
et al., in preparation). The fraction of the young 
stellar mass over the total mass is well constrained, 
with values fySP < 0.26 for all the galaxies, except 
for NGC1705, NGC2798 and NGC4631, for which 
fySP is unconstrained. The mean value for the pa- 
rameter a of the Dale & Helou (2002) templates 
is a = 2.44 (in agreement with Noh et al. (2009); 
we note that, while low values of a are well con- 
strained, no constraints are found on high values, 
with a — 3 compatible at 68% c.l. for all galaxies 
expect NGC2798, for which a < 1.77 at 68% c.l.; 
weak constraints on high values for a are due to the 
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degeneracy of the dust emission models for wave- 
lengths shortwards of the emission peak, so that, 
for a > 2.5; the flux ratio fi/(60)/iTO/f^(100)/im is 
almost constant (see Noll et al. (2009). 

• Good constraints can be derived for the mass 
dependent parameters Mgtar, SFR, Lboi and 
Ldust as shown in Table 3. The Pearson cor- 
relation coefficient with results from Noll et 
al. (2009) is r = 0.99, 0.98, 0.94, 1.0, 1.0, 0.94 
for Mstar, SFR, tD4000, Lboi, Ldust, Afuv respec- 
tively. In Figure 8 we show a quantitative com- 
parison with the analysis performed by Noll et al. 
(2009) with the original code CIGALE by plotting 
the ratio, Q = (CIGALEMC - CIGALE)/ctcigALEMC> 
where CIGALEMC and CIGALE refer to the mean 
value of the parameters quoted in Table 3 for this 
work and Noll et al. (2009) respectively. Possible 
systematic differences between the results of both 
methods can be studied by considering mean, stan- 
dard deviation and skewness for the parameters of 
interest. As we can see from Table 4, no signifi- 
cant difference between results in this paper and 
in Noll et al. (2009) is found for SFR, log(tD4ooo), 
log(Lboi), and Apuv, with each value compatible 
with at less or about 68% c.l.; however, there 
is some skewness associated to SFR, tD4ooo and 
Lboi- In case of Mgtar and Ldust, the mean val- 
ues are lower from at about 95% c.l., indicating 
the existence of some offset between results in this 
paper and in Noll et al. (2009). A further compar- 
ison of the uncertainties between our results and 
Noll et al. (2009) shows that we are able to put 
stronger constraints on the galaxy SFR, while our 
constraints on Mgtar, Lboi, Ldust, Apuv are weaker. 
Differences between our results and results in Noll 
et al. (2009) can be due to the different choice of 
the input parameter space in terms of both num- 
ber and type of parameters used in the analysis; in 
particular, we consider a wider range of variation 
for a, 6, Aysp, t2 and fAGN, while we decide to fix 
the parameter t2 (see Table 3 of Noll et al. (2009) 
for their choice of the parameter space); however, 
it is reassuring to see how different theoretical as- 
sumptions lead to compatible results. 

• A frequency dependent analysis shows which 
bands mostly contribute to the total for each 
galaxy. We introduce the averaged x^- 

J i=l "i 

here Th; and d] are respectively the theoretical best 
fit SEDs and the data points for the i-th galaxies at 
the j-th frequency and the sum runs over our sam- 
ple of galaxies. From Figure 7 we see that the code 
is able to find a good agreement with the data for 
all the bands considered. The frequency dependent 
is particularly low for bands with large relative 
errors as UV GALEX bands and optical B, V, R, I 
bands. The large uncertainty found for the MIPS 
70/im filter is mainly due to the galaxies NGC1715 
and NGC5866; the code is not able to find a good 
fit for the 70/im filter and the reduced chi-squared 



-2 



-4 



-s 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
^AGN 

Fig . 6 . — Two-dimensional marginalized distribution showing the 
68% and %95 c.l. contours for log(SFR) vs Jagn for NGC1404; 
the double peaked likelihood is clearly visible and fAGN is mostly 
unconstrained. 

is high {x'^/d.o.f = 7.8 for both galaxies, see Table 
3). A previous analysis by Cannon et al. (2006), 
based on Spitzer observations, found similar results 
for NGC1705, showing in particular that the mod- 
els of Li & Draine (2001, 2002) give a better fit to 
IR data than the Dale & Helou (2002) models used 
in this paper (see Figure 3 of Cannon et al. (2006)). 
Interestingly enough, NGC1715 and NGC5866 are, 
respectively, the only irregular and SO galaxies in 
our sample. 

In general, we find a small trend of worse fitting for 
galaxies with small SFR values. Galaxies with the 
lowest values of the do not show peculiar proper- 
ties, with total stellar masses of about IO^^Mq and 
SFRs between 0.1 and IOMq /yr, typical of nearby 
spiral galaxies; an exception is the dwarf galaxy 
NGC4625 (reduced x^ld.o.f = 1.2), characterized 
by smaller values for both the stellar mass and lu- 
minosity. Three other dwarf galaxies (NGC1705, 
NGC2976 and NGC5474) are clearly identified in 
our sample by looking at their smaller values for 
both the luminosity and the stellar masses respect 
to the rest of the sample. Finally, in order to check 
for a possible correlation between the filter set used 
and the goodness of the fit, we also calculated the 
mean x^ for the subset of galaxies for which SDSS 
filters are available respect to the subset of galax- 
ies with optical Dale et al. filters; we found no 
correlation, with x ~ 2.5 in both cases. 

6. CONCLUSIONS 

In this paper we have introduced a MCMC sampling 
method for the astrophysical parameter estimation from 
SED fitting with CIGALE. We have shown the following 
advantages of our modified CIGALEMC code over the usual 
grid-based CIGALE: 

• its efficiency, in terms of CPU time, through the 
Metropolis-Hastings algorithm. Most of the sam- 
pling points are drawn in the region where the pos- 
terior probability is high, while in the grid-based 
approach all regions are sampled in the same way. 
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TABLE 3 

Mean values and 68% c.l. marginalized results for some parameters related to the SINGS sample considered. 



ID 


Type 


log Mstar 


log SFR 


log tD4000 


log ibol 


log idust 


^FUV 


X^/d.o.f. 






[M©] 


[Me/yr] 


[Gyr] 


[L0] 


[L0] 


[mag] 







SAc 


9 


52 


-|- 


Q 


10 


-0.78 


-|- 


Q 


XI 


-0 


14 


-|- 


Q 


18 


9.46 


-|- 


Q 


05 


8.62 


-|- 


Q 


XX 


0.45 


-|- 


Q 


16 


2 


5 


NOP n^84 


E4 


1^ 


41 


-|- 


Q 


04 


-1.3 


-|- 


Q 


25 


Q 


99 


-|- 


Q 


02 


11.09 


-|- 


Q 


18 


10.22 


-|- 


X 


14 


4.47 


-|- 


4 


41 


2 




NOP nq^'i 


SABd 


9 


93 


-|- 


Q 


17 


0.18 


-|- 


Q 


12 


-0 


43 


-|- 


Q 


10 


10.24 


-|- 


Q 


07 


9.64 


-|- 


Q 


09 


0.6 


-|- 


Q 


21 


2 




NOC inQ7 


SBb 




19 


-|- 





XI 


0.9 


-|- 





XX 


-0 


15 


-|- 





17 


11.15 


-|- 





04 


10.77 


-|- 





08 


1.85 


-|- 





51 





5 


NOP 1 9Q1 


SBa 




14 


-|- 


Q 


05 


-0.72 


-|- 


Q 


23 


Q 


86 


-|- 


Q 


08 


10.68 


-|- 


Q 


04 


9.24 


-|- 


Q 


16 


0.98 


-|- 


Q 


66 


X 


9 


NOP 1 '^1 fi 


SABO 


12 


01 


-|- 


Q 


05 


-0.06 


-|- 


Q 


28 


Q 


91 


-|- 


Q 


07 


11.54 


-|- 


Q 


04 


10.07 


-|- 


Q 


12 


1.49 


-|- 


X 


07 


4 


5 


NPP 1 404 


El 




52 


-|- 


Q 


04 


-1.0 


-|- 


Q 


20 


Q 


98 


-|- 


Q 


02 


11.25 


-|- 


Q 


21 


10.76 


-|- 


Q 


58 


5.58 


-|- 


5 


07 


X 


4 


NOP 1 ^12 


SBab 


10 


34 


-|- 





09 


-0.28 


-|- 





14 





17 


-|- 





32 


10.13 


-|- 





04 


9.39 


-|- 





09 


0.86 


-|- 





30 


X 


8 






10 


88 


-|- 


Q 


XI 


0.9 


-|- 


Q 


10 


-0 


34 


-|- 


Q 


08 


11.02 


-|- 


Q 


05 


10.58 


-|- 


Q 


08 


1.13 


-|- 


Q 


34 


X 




NOP 1 70'^ 


Am 




20 


-|- 


Q 


20 


-1.15 


-|- 


Q 


15 


-0 


62 


-|- 


Q 


25 


8.81 


-|- 


Q 


10 


7.67 


-|- 


Q 


12 


0.11 


-|- 


Q 


05 






NHP 97QR 


SBa 


10 


03 


-|- 


Q 


18 


0.61 


-|- 


Q 


07 


-0 


54 


-|- 


Q 


06 


10.6 


-|- 


Q 


05 


10.47 


-|- 


Q 


06 


4.48 


-|- 


Q 


64 


2 


Q 


NCIC 2841 


SAb 


10 


92 


-|- 





06 


-0.3 


-|- 





17 





62 


-|- 





10 


10.54 


-|- 





03 


9.57 


-|- 





XX 


1.04 


-|- 





46 


X 


3 


NPP 9q7fi 


SAc 


9 


33 


-|- 


Q 


09 


-0.8 


-|- 


Q 


08 


-0 


28 


-|- 


Q 


06 


9.37 


-|- 


Q 


04 


8.87 


-|- 


Q 


09 


1.11 


-|- 


Q 


28 


3 


4 


NPP '^n^i 


SAab 


10 


96 


-|- 


Q 


06 


-0.16 


-|- 


Q 


14 


Q 


56 


-|- 


Q 


10 


10.59 


-|- 


Q 


03 


9.51 


-|- 


Q 


XX 


0.63 


-|- 


Q 


27 


X 


g 


NOP S1S4 




10 


14 


-|- 





08 


0.11 


-|- 





07 


-0 


33 


-|- 





07 


10.24 


-|- 





03 


9.68 


-|- 





09 


0.78 


-|- 





23 


4 


4 


NOP siqn 


SAap 


10 


87 


-|- 





04 


-0.65 


-|- 





35 





75 


-|- 





15 


10.44 


-|- 





03 


9.62 


-|- 





10 


3.16 


-|- 


X 


23 


X 


4 


NPP '^1 QS 


SBc 


10 


01 


-|- 


Q 


08 


-0.02 


-|- 


Q 


07 


-0 


33 


-|- 


Q 


07 


10.12 


-|- 


Q 


04 


9.55 


-|- 


Q 


09 


0.73 


-|- 


Q 


21 


2 


2 


NPP S'^'il 


SBb 


10 


58 


-|- 


Q 


08 


-0.02 


-|- 


Q 


12 


Q 


14 


-|- 


Q 


23 


10.38 


-|- 


Q 


04 


9.82 


-|- 


Q 


09 


1.36 


-|- 


Q 


41 


X 


X 


NPP '^^91 




10 


95 


-|- 


Q 


08 


0.41 


-|- 


Q 


14 


Q 


06 


-|- 


Q 


23 


10.76 


-|- 


Q 


04 


10.28 


-|- 


Q 


10 


2.06 


-|- 


Q 


51 


Q 


9 


NOP Sfi21 


SAd 


10 


04 


-|- 





12 


0.14 


-|- 





10 


-0 


38 


-|- 





08 


10.24 


-|- 





05 


9.82 


-|- 





09 


1.15 


-|- 





35 


X 


4 


NPP '^627 


SABb 


10 


80 


-|- 


Q 


10 


0.55 


-|- 


Q 


XX 


-0 


21 


-|- 


Q 


08 


10.77 


-|- 


Q 


04 


10.38 


-|- 


Q 


09 


2.0 


-|- 


Q 


40 


X 




NPP 4^'^fi 


SARbr 


10 


89 


-|- 


Q 


12 


1.0 


-|- 


Q 


10 


-0 


39 


-|- 


Q 


08 


11.11 


-|- 


Q 


04 


10.79 


-|- 


Q 


08 


1.65 


-|- 


Q 


40 


X 


2 


NPP 4^^q 




10 




-|- 


Q 


16 


0.43 


-|- 


Q 


08 


-0 


47 


-|- 


Q 


06 


10.46 


-|- 


Q 


04 


9.84 


-|- 


Q 


09 


0.52 


-|- 


Q 


15 


3 


3 


NPP 4^fiq 


Q A Rah 




38 


-|- 


Q 


05 


0.33 


-|- 


Q 


17 


Q 


53 


-|- 


Q 


12 


11.03 


-|- 


Q 


03 


10.26 


-|- 


Q 


10 


1.66 


-|- 


Q 


56 


4 


Q 


NPP 4^7q 


SABb 


11 


42 


-|- 


Q 


06 


0.18 


-|- 


Q 


21 


Q 


63 


-|- 


Q 


XX 


11.03 


-|- 


Q 


03 


10.13 


-|- 


Q 


XX 


1.5 


-|- 


Q 


60 


X 


3 


NPP 4^q4 


SAa 


11 


73 


-|- 


Q 


05 


-0.44 


-|- 


Q 


27 


Q 


94 


-|- 


Q 


06 


11.26 


-|- 


Q 


04 


9.7 


-|- 


Q 


12 


1.29 


-|- 


Q 


97 


X 


4 


NGG 4625 


SABmp 


9 


18 


± 





10 


-0.77 


± 





08 


-0 


37 


± 





08 


9.33 


± 





04 


8.79 


± 





10 


0.75 


± 





22 


1 


2 


NOG 4631 


SBd 


10 


08 


± 





17 


0.76 


± 





07 


-0 


59 


± 





06 


10.73 


± 





04 


10.45 


± 





08 


1.39 


± 





36 


2 


9 


NGG 4725 


SABab 


11 


34 


± 





07 


0.31 


± 





15 





51 


± 





11 


11.0 


± 





03 


10.01 


± 





12 


0.7 


± 





29 


2 


3 


NGG 4736 


SAab 


10 


75 


± 





07 


0.08 







12 





21 


± 





30 


10.51 


± 





03 


9.84 







10 


1.2 


± 





36 


1 


7 


NGG 4826 


SAab 


10 


77 


± 





05 


-0.44 


± 





21 





62 


± 





12 


10.38 


± 





03 


9.54 


± 





11 


1.78 


± 





68 


2 


6 


NGG 5033 


SAc 


10 


68 







10 


0.43 


± 





11 


-0 


20 







09 


10.65 







04 


10.25 


± 





09 


1.7 


± 





42 


1 


2 


NGG 5055 


SAbc 


10 


88 


± 





08 


0.45 


± 





11 


-0 


06 


± 





14 


10.75 


± 





04 


10.27 


± 





09 


1.7 


± 





41 


2 


2 


NGG 5194 


SABbc 


10 


74 


± 





11 


0.85 


± 





09 


-0 


39 


± 





08 


10.94 


± 





04 


10.55 


± 





09 


1.3 


± 





33 





8 


NGG 5195 


SBOp 


10 


76 


± 





04 


-0.42 


± 





27 





61 


± 





15 


10.38 


± 





04 


9.61 


± 





15 


2.55 


± 





92 


5 


6 


NGG 5474 


SAcd 


9 


22 


± 





14 


-0.53 


± 





10 


-0 


43 


± 





07 


9.52 


± 





04 


8.52 


± 





12 


0.19 


± 





07 


5 


6 


NGG 5713 


SABbcp 


10 


54 


± 





19 


0.83 


± 





11 


-0 


44 


± 





09 


10.89 


± 





05 


10.68 


± 





08 


2.74 


± 





56 


1 


1 


NGG 5866 


SO 


10 


88 


± 





04 


-1.01 


± 





36 





87 


± 





10 


10.42 


± 





03 


9.2 


± 





14 


2.56 


± 


1 


26 


7 


8 


NGG 7331 


SAb 


11 


31 


± 





11 


0.78 







15 





08 







32 


11.13 







04 


10.73 


± 





09 


2.62 


± 





64 


1 


2 



10.00 



1.00 



0.10 



0.01 
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Parameters 


mean 


VVar 


Skewness 


logio Mstar 


-1.0 


0.52 


-0.1 


logio SFR 


0.29 


0.89 


1.79 


logtD4000 


1.2 


2.8 


2.6 


logio Lbol 


-0.45 


0.36 


1.17 


logio Ldust 


-0.58 


0.25 


0.01 


Apuv 


-0.69 


0.62 
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Fig. 7. — The averaged frequency dependent for our SINGS 
sample shows a general good fit for all the bands used; the large 
dispersion for the TOfim MIPS filter is mainly due to NGG1705 and 
NGG5866. Error bars are calculated as standard deviations. 



TABLE 4 

Values of the first three moments of the 

(CIGALEMC — CIGALE)/aQlQf^l^^f^Q distributions FOR SOME 
parameters OF INTEREST. 



Moreover, marginalized one-dimensional probabil- 
ity distributions for the parameters of interest are 
calculated by simply counting the number of sam- 
ples within a binned range of parameter values, 
the density of sampling points being proportional 
to the posterior probability. It is hard to do the 
same with the usual grid approach, since the in- 
tegration calculation scales exponentially with the 
number of dimensions. The analysis of a mock 
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Fig. 8. — Ratio Q = (CIGALEMC — CIGALE)/o-cigalemc between the mean values estimated in this work and in Noll et al. (2009) 
for some parameters quoted in Table 3. 



sample shows that CIGALEMC needs 20 times less 
points than CIGALE to reach convergence for a 
given galaxy but, in general, results depend on both 
the number of parameters and the prior used. 

• its accuracy; degeneracies between parameters are 
easily found, convergence criteria (already imple- 
mented in the code) ensure that statistical quanti- 
ties for the parameters of interest are robustly de- 
termined and cross-checks through statistical anal- 
ysis of mock catalogs are not necessary. 

• its "user friendly" characteristics. The user does 
not need to decide a priori the number density 
of samples for each region, trying to find a com- 
promise between the accuracy of the results and 



the speed of the code: only the prior range must 
be chosen in advance for CIGALEMC; in fact the 
Metropolis-Hastings algorithm automatically sam- 
ples adequately the posterior probability according 
to its values. 

Our code will be available ver y soon at this web address: 



http: / /www. oamp.fr/cigale/ 
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